SETUP

Install and load packages

knitr::opts_chunk$set(echo = TRUE)

# install.packages(
#       "ape", "phyloregion", "canaper",            # phylo libraries
#       "terra", "geosphere", "prioritizr", "gdm",  # spatial libraries
#       "pals", "vegan", "tidyverse", "Rsymphony")  # other utilities

# install.packages("canaper", repos = "https://ropensci.r-universe.dev")

library(ape)
library(tidyverse)
library(terra)
library(phyloregion)
library(canaper)
library(prioritizr)
library(pals)
library(vegan)
library(geosphere)

select <- dplyr::select

Define a few helper functions for later use

# function to clean a spatial phylogenetic dataset
clean_dataset <- function(x){

      # remove spaces from taxon names
      x$tree$tip.label <- str_replace(x$tree$tip.label, " ", "_")
      colnames(x$comm) <- str_replace(colnames(x$comm), " ", "_")

      # convert non-finite values to formal absences
      x$comm[!is.finite(x$comm)] <- 0

      # remove taxa with no occurrences
      ghost <- colSums(x$comm, na.rm = T) == 0
      if(any(ghost)){
            message("removing ", sum(ghost), " taxa with no occurrences")
            x$comm <- x$comm[, !ghost]
      }

      # remove non-overlapping taxa
      taxa <- intersect(x$tree$tip.label, colnames(x$comm)[colSums(x$comm, na.rm = T) > 0])
      if(length(taxa) == 0) stop("no overlapping taxon names")
      prune <- setdiff(x$tree$tip.label, taxa)
      if(length(prune > 0)){
            message("removing ", length(prune), " tips from tree")
            x$tree <- drop.tip(x$tree, prune)
      }
      drop <- setdiff(colnames(x$comm), taxa)
      if(length(drop > 0)){
            message("removing ", length(drop), " taxa from community matrix")
            x$comm <- x$comm[, taxa]
      }

      # remove sites with no taxa
      occupied <- rowSums(x$comm) > 0
      if(sum(occupied > 0)) message("removing ", sum(occupied), " sites with no occurrences")
      x$comm <- x$comm[occupied,]
      x$xy <- x$xy[occupied,]

      # name rows
      rownames(x$comm) <- paste0("cell_", 1:nrow(x$comm))

      return(x)
}

# a set of functions to construct a community matrix
# with a column for every branch of a phylogeny, not just every tip
parentProb <- function(x) 1 - prod(1 - x)
build_clade_range <- function(e, phylo, comm){
      node <- phylo$edge[e,2]
      if(node <= length(phylo$tip.label)){
            otu <- phylo$tip.label[node]
            prob <- comm[,otu]
      } else{
            clade <- extract.clade(phylo, node)
            otu <- clade$tip.label
            prob <- apply(comm[,otu], 1, parentProb)
      }
      return(prob)
}
build_clade_ranges <- function(tree, comm){
      sapply(1:nrow(tree$edge),
             build_clade_range, phylo=tree, comm=comm)
}

# helpher function to plot a phylogeny with edge colors
phyplot <- function(tree,
                   connect = NULL, hl = "orange", bg = "black",
                   value = NULL, col = c("black", "blue", "red", "orange"),
                   ...){

      clr <- rep(bg, length(tree$edge.length))
      if(!is.null(connect)){
            clr[which.edge(tree, connect)] <- hl
      }

      if(!is.null(value)){
            n <- min(c(20, length(tree$edge.length)))
            pal <- colorRampPalette(col)(n)
            clr <- pal[cut(value, n)]
            if(sd(value[is.finite(value)]) == 0) clr <- "black"
      }

      plot(tree, edge.color = clr, ...)
}

# helper function to ggplot a raster map
carto <- function(d, v){
      p <- ggplot() +
            geom_tile(data = d,
                      aes_string("x", "y", fill = v)) +
            xlim(range(d$x) + diff(range(d$x))/20 * c(-1, 1)) +
            ylim(range(d$y) + diff(range(d$y))/20 * c(-1, 1)) +
            coord_fixed(ratio = 1.2) +
            theme_void() +
            theme(legend.position = "top")
      if(inherits(d[[v]], "numeric")) p <- p + scale_fill_viridis_c()
      p
}

WARMUP

Phylogenies

Let’s start with a quick primer on phylogenetic trees in R. We’ll use ape, the core R phylogenetics library. (Other libraries to be aware of for working with phylogenies include phytools, picante, and ggtree, among others.)

# load a tree and plot it
tree <- read.tree("../../data/mishler_2014/mishler_2014_acacia.tre")
plot(tree, type = "fan", show.tip.label = F)

# explore the data structure of a phylogeny
# class(tree)
# print(tree)
# str(tree)
# ?plot.phylo
# plot(tree, show.tip.label = F, type = "radial")

# plot the branches connecting a set of tips (e.g. a local community) --
# PD is the lengths of these branches (plus the edges connecting them to the root)
community <- sample(tree$tip.label, 5)
phyplot(tree, connect = community, 
        show.tip.label = F)

Exercises:

  • Load and plot the mammal phylogeny (Upham et al. 2019).
  • Plot the subtree representing the marsupials. The most recent common ancestor of the two species listed below is the crown node of all marsupials; you can find it using getMRCA() and prune it from the full tree using extract.clade().
  • Plot the full mammal tree with marsupials highlighted in red. You can use phytools::getDescendants() to index all the descendants of the marsupial MRCA.
# Load and plot the mammal phylogeny
marsupials <- c("Lestoros_inca", "Lagostrophus_fasciatus")

# Plot the subtree representing the marsupials
tree <- read.tree("../../data/upham_2019/upham_2019_tree0000.tre")
mrca <- getMRCA(tree, marsupials)
extract.clade(tree, mrca)
## 
## Phylogenetic tree with 362 tips and 361 internal nodes.
## 
## Tip labels:
##   Rhyncholestes_raphanurus, Lestoros_inca, Caenolestes_fuliginosus, Caenolestes_sangay, Caenolestes_caniventer, Caenolestes_condorensis, ...
## 
## Rooted; includes branch lengths.
# Plot the full mammal tree with marsupials highlighted in red
desc <- phytools::getDescendants(tree, mrca)
clr <- rep("black", nrow(tree$edge))
clr[which.edge(tree, desc)] <- "red"
plot(tree, edge.color = clr, show.tip.label = F, type = "fan")

Spatial data in R

For a spatial phylogenetic analysis, we need to define geographic areas of interest, within which organisms are considered part of the same community. These areas can be grid cells or can be irregular polygons, but we’ll focus on the former here as the latter are conceptually similar but more complex to work with.

We’ll be working with raster data and point data. The raster data comprise a spatial grid of presences and absences. Point data represent occurrence localities, which we’ll convert to grids of presences and absences before analysis.

Here’s a quick demo of loading and visualizing some spatial data.

# raster data, read in with the terra package
r <- rast("../../data/jetz_2012/bbs_occ.tif")
plot(r[[1:4]], col = "green")

plot(sum(r, na.rm = T))

# point data, read in as a data frame
p <- read_csv("../../data/mishler_2014/mishler_2014_acacia_points.csv")
## Rows: 132248 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): Species, taxa
## dbl (4): Longitude, Latitude, x_metres_EPSG_3577_Albers_Equal_Area, y_metres...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# glimpse(p)
ggplot(p, aes(Longitude, Latitude)) + 
      geom_point()

# we can "rasterize" points by rounding the coordinates
# (there are fancier approaches but we'll keep it simple here)
pr <- p %>% mutate(x = round(Longitude), y = round(Latitude))
ggplot(pr, aes(x, y)) + 
      geom_raster()

# taking this a step farther, lets' convert our rounded coordinates 
# for one species into to a proper raster object
rp <- pr %>%
      filter(taxa == taxa[1000]) %>%
      dplyr::select(x, y) %>%
      distinct() %>%
      mutate(present = T) %>%
      rast()
plot(rp, col = "black")

Exercises:

  • Load the terrestrial mammal spatial data using rast()
  • Plot a map of the range of the gray wolf (Canis lupus).
  • Plot a map of canid species richness (all species in the genus Canis).
  • Plot a map of overall mammal species richness.
r <- rast("../../data/upham_2019/iucn_mammals.tif")
r %>% subset("Canis lupus") %>% plot()

r %>% subset(grep("Canis ", names(r))) %>% sum() %>% plot()

r %>% sum() %>% plot()

SPATIAL PHYLOGENETIC DATA

At a minimum, every spatial phylogenetic dataset consists of a phylogeny and a spatial dataset representing occurrence locations of the terminal taxa. We’ll use four different datasets as examples.

To keep things tidy, let’s start by loading the data for each one, and getting the pieces into a standardized format. We’ll package each dataset in a list with three pieces:

To avoid issues with some of the algorithms used later on, we need to do a few things to clean up the datasets, including removing taxa that aren’t present in both the spatial and phylogenetic datasets, and removing unoccupied sites from the community matrix. To do this we’ll use a function called clean_dataset() defined in functions.R.

Cailfornia vascular plants (from Thornhill et al. 2017):

# phylogeny
tree <- read.nexus("../../data/thornhill_2017/thornhill_2017_chronogram.nex")

# spatial data
occ <- read_csv("../../data/thornhill_2017/thornhill_2017_occ.csv")
## Rows: 1394170 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): otu
## dbl (2): x, y
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
glimpse(occ)
## Rows: 1,394,170
## Columns: 3
## $ x   <dbl> -121.8009, -121.0989, -122.2753, -121.6883, -121.2134, -119.5225, …
## $ y   <dbl> 37.39750, 40.07583, 40.74150, 37.21714, 40.02400, 37.44472, 33.533…
## $ otu <chr> "Delphinium", "Juncus", "Cyperus__Lipocarpha", "Cryptantha", "Cype…
# # compare species lists
all.equal(sort(tree$tip.label),
          sort(unique(occ$otu)))
## [1] TRUE
# summarize to lat-long grid
rnd <- function(x, y) round(x/y)*y # round x to the nearest y
occ <- occ %>%
      mutate(x = rnd(x, .5),
             y = rnd(y, .5)) %>%
      distinct()

# construct community matrix
comm <- occ %>%
      mutate(pres = 1) %>%
      spread(otu, pres, fill = 0)
xy <- select(comm, x, y)
comm <- select(comm, -x, -y) %>%
      as.matrix()

# package data
flora <- list(tree = tree,
              comm = comm,
              xy = xy) %>%
      clean_dataset()
## removing 221 sites with no occurrences

Australian acacia (from Mishler et al. 2014):

# phylogeny
tree <- read.tree("../../data/mishler_2014/mishler_2014_acacia.tre")

# spatial data
occ <- read_csv("../../data/mishler_2014/mishler_2014_acacia_points.csv")
## Rows: 132248 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): Species, taxa
## dbl (4): Longitude, Latitude, x_metres_EPSG_3577_Albers_Equal_Area, y_metres...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# construct community matrix
occ <- occ %>%
      select(otu = taxa, x = Longitude, y = Latitude) %>%
      mutate(x = rnd(x, 1), # summarize to 0.5 degree lat-lon grid
             y = rnd(y, 1)) %>%
      distinct()
comm <- occ %>%
      mutate(pres = 1) %>%
      spread(otu, pres, fill = 0)
xy <- select(comm, x, y)
comm <- select(comm, -x, -y) %>%
      as.matrix()

# package data
acacia <- list(tree = tree,
               comm = comm,
               xy = xy) %>%
      clean_dataset()
## removing 797 sites with no occurrences

American birds (phylogeny from Jetz et al. 2012, spatial data from US Breeding Bird Survey):

# phylogeny (one of 10,000 posterior trees from Jetz 2012)
tree <- read.nexus("../../data/jetz_2012/jetz_2012_tree0001.nex")

# spatial data
r <- rast("../../data/jetz_2012/bbs_occ.tif")

# construct community matrix
comm <- values(r)
dim(comm)
## [1] 377 442
image(t(comm))

# spatial coordinates
xy <- crds(r[[1]], df = T, na.rm = F)

# package data
birds <- list(tree = tree,
              comm = comm,
              xy = xy) %>%
      clean_dataset()
## removing 242 sites with no occurrences

Mammals (phylogeny from Upham et al. 2019, spatial data from IUCN):

# phylogeny (one of 10,000 posterior trees from Upham 2019)
tree <- read.tree("../../data/upham_2019/upham_2019_tree0000.tre")

# spatial data
r <- rast("../../data/upham_2019/iucn_mammals.tif")

# construct community matrix
comm <- values(r)

# spatial coordinates
xy <- crds(r[[1]], df = T, na.rm = F)

# package data
mammals <- list(tree = tree,
                comm = comm,
                xy = xy) %>%
      clean_dataset()
## removing 734 tips from tree
## removing 446 taxa from community matrix
## removing 2598 sites with no occurrences

Exercises:

  • Plot the phylogeny for each of the datasets.
  • Plot the spatial coordinates for each of the datasets.
  • Identify the number of taxa and number of sites in each dataset.
mammals$tree %>% plot(show.tip.label = F)

mammals$xy %>% plot()

dim(mammals$comm)
## [1] 2598 5177

ALPHA DIVERSITY

Alpha phylogenetic diversity (and related metrics like phylogenetic endemism) represent the amount of evolutionary history present in a single community. We’ll look at functions from some existing packages that can calculate phylodiversity metrics, and reproduce these with some of our own calculations. We’ll also cover methods for significance testing, and explore the sensitivity of these tests to spatial scale. And we’ll take a look at how alternative evolutionary branch lengths influence results.

Pre-packaged functions

ds <- birds

# with phyloregion library
scomm <- ds$comm %>% Matrix(sparse = TRUE)
# PD(scomm, ds$tree)
pr <- ds$xy %>%
      mutate(pd = PD(scomm, ds$tree),
             pe = phylo_endemism(scomm, ds$tree),
             we = weighted_endemism(scomm))
pr %>% ggplot(aes(x, y, fill = pd)) + geom_raster()

pr %>% carto("pd")

pr %>% carto("pe")

# with canaper library
# ?cpr_rand_test
cpr <- cpr_rand_test(ds$comm, ds$tree, 
                     null_model = "curveball",
                     n_reps = 1) %>%
      bind_cols(ds$xy)
# glimpse(cpr)
# names(cpr)
cpr %>% carto("pd_obs")

cpr %>% carto("pe_obs")

# plot a local community on the phylogeny
par(mfrow = c(1, 2))
for(fun in list(min, max)){
      cell <- cpr %>% filter(pd_obs == fun(pd_obs)) %>% rownames()
      spp <- colnames(ds$comm)[ds$comm[cell,] == 1]
      phyplot(ds$tree, connect = spp, 
              type = "fan", show.tip.label = F,
              main = paste("taxa found in", cell))
}

Manual approach

Canaper and phyloregion are pretty great, but how are these spatial biodiversity metrics actually being calculated under the hood? Let’s reproduce some of these results with a slightly more manual approach.

# build community matrix including not just tips as above,
# but also every clade at every level
pcomm <- build_clade_ranges(ds$tree, ds$comm) # from functions.r
# dim(ds$comm)
# dim(pcomm)
image(t(pcomm))

# clade-level metrics
range_size <- colSums(pcomm)
branch_length <- ds$tree$edge.length / sum(ds$tree$edge.length)
phyplot(ds$tree, value = branch_length,
        show.tip.label = F, type = "fan", main = "branch length")

phyplot(ds$tree, value = log(1 / range_size),
        show.tip.label = F, type = "fan", main = "log endemism")

# PD and PE are just weighted summaries across this pcomm matrix:
pd <- pcomm %>%
      apply(1, function(x) x * branch_length) %>%
      colSums()
ds$xy %>% mutate(pd = pd) %>% carto("pd")

pe <- pcomm %>%
      apply(1, function(x) x * branch_length / range_size) %>%
      colSums()
ds$xy %>% mutate(pe = pe) %>% carto("pe")

Exercises:

  • Compute these metrics with one of the other datasets.
  • Species richness can be considered PD on a star phylogeny. Prove it.
  • Manual calculations make it clear that these diversity metrics are simply weighted summaries of the pcomm matrix, which invites any number of additional weighting schemes. Make a version weighted by: taxon endangerment; site endangerment; or taxon-site habitat importance (using fake simulated data for these variables).
# SR == PD on star phylogeny
ds <- birds
star <- ds$tree
tip <- ! star$edge[,2] %in% star$edge[,1]
star$edge.length <- rep(0, length(star$edge.length))
star$edge.length[tip] <- 1
plot(star, type = "fan", show.tip.label = F)

sr <- rowSums(ds$comm)
pd <- PD(scomm, star)
plot(sr, pd)

# taxon & site endangerment
# (this would represent locations of treatened sites with high concentrations of range-restricted, endangered taxa)
endangerment <- runif(ncol(pcomm)) # fake random data
site_threat <- runif(nrow(pcomm)) # fake random data
tpe <- pcomm %>%
      apply(1, function(x) x * branch_length / range_size * endangerment) %>%
      apply(1, function(x) x * site_threat) %>%
      rowSums()
ds$xy %>% mutate(tpe = tpe) %>% carto("tpe")

Significance testing

To get a sense for how noteworthy an observed high or low diversity metric is, it can be helpful to compare the observed value to a null distribution. Let’s demonstrate this using the acacia dataset. The cpr_rand_test() function randomly permutes the community matrix many times, calculating the diversity of each location each time to derive a distribution of null diversity values for each location.

CANAPE (Mishler et al. 2014) is a special application of this randomization approach, used to identify areas of neo-endemism and paleo-endemism.

ds <- acacia

# note that this is a phylogram, not a chronogram
# ds$tree %>% plot(show.tip.label = F)

# phylodiversity randomizations
cpr <- cpr_rand_test(ds$comm, ds$tree, 
                     null_model = "curveball",
                     n_reps = 101, n_iterations = 10000) %>%
      bind_cols(ds$xy)

# raw PD patterns here are very different from randomized quantiles
cpr %>% carto("pd_obs")

cpr %>% carto("pd_obs_p_upper")

cpr %>%
      ggplot(aes(pd_obs, pd_obs_p_upper)) +
      geom_point()

cpr %>% carto("rpd_obs_p_upper")

## CANAPE ##

# classify canape categories
end <- cpr %>% cpr_classify_endem() 

# map
end %>%
      carto("endem_type") +
      scale_fill_manual(
            values = c("mediumpurple1", "red", "gray80", "blue", "darkorchid4"),
            breaks = c("mixed", "neo", "not significant", "paleo", "super"))

# scatterplot of canape components
end %>%
      ggplot(aes(pe_alt_obs, pe_obs,
                 color = endem_type)) +
      geom_point() +
      scale_color_manual(
            values = c("mediumpurple1", "red", "gray80", "blue", "darkorchid4"),
            breaks = c("mixed", "neo", "not significant", "paleo", "super"))

Null models can be constructed in many different ways, and these differences can have a major effect on results. We won’t get too far into the weeds here, but let’s illustrate how a few different null models affect PD significance values for this dataset:

# ?vegan::commsim
curveball <- cpr %>% 
      mutate(model = "curveball")
r0 <- cpr_rand_test(ds$comm, ds$tree, null_model = "r0", metrics = "pd", 
                    n_reps = 101) %>%
      bind_cols(ds$xy) %>%
      mutate(model = "r0")
c0 <- cpr_rand_test(ds$comm, ds$tree, null_model = "c0", metrics = "pd", 
                    n_reps = 101) %>%
      bind_cols(ds$xy) %>%
      mutate(model = "c0")

bind_rows(curveball, r0, c0) %>%
      mutate(pd_sig = case_when(pd_obs_p_upper > .95 ~ "high",
                                pd_obs_p_upper < .05 ~ "low",
                                TRUE ~ "no")) %>%
      carto("pd_sig") + 
      facet_wrap(~model, nrow = 1) +
      scale_fill_manual(values = c("red", "blue", "gray"))

Before moving on from randomizations, it’s important to note that the size of the geographic region in an analysis can strongly influence the results, because it alters the species pool to which a site is compared during randomization. Let’s look at three different extents:

# a function to subset the sites in a dataset
filter_sites <- function(ds, selection){
      dsf <- ds
      dsf$comm <- dsf$comm[selection, ]
      dsf$xy <- dsf$xy[selection, ]
      dsf<- clean_dataset(dsf)
      return(dsf)
}

# filter mammals data to three nested regions of different sizes
americas <- mammals %>% 
      filter_sites(mammals$xy$x < -34)
## removing 3415 taxa with no occurrences
## removing 3415 tips from tree
## removing 914 sites with no occurrences
s_america <- mammals %>% 
      filter_sites(between(mammals$xy$x, -82, -34) 
                   & mammals$xy$y < 13)
## removing 4027 taxa with no occurrences
## removing 4027 tips from tree
## removing 227 sites with no occurrences
amazon <- mammals %>% 
      filter_sites(between(mammals$xy$x, -75, -50) 
                   & between(mammals$xy$y, -15, 5))
## removing 4562 taxa with no occurrences
## removing 4562 tips from tree
## removing 56 sites with no occurrences
# compute pd significance for each dataset
cpr_americas <- cpr_rand_test(
      americas$comm, americas$tree, 
      null_model = "curveball", metrics = "pd", n_reps = 21) %>%
      bind_cols(americas$xy) %>%
      mutate(region = "americas")
cpr_s_america <- cpr_rand_test(
      s_america$comm, s_america$tree, 
      null_model = "curveball", metrics = "pd", n_reps = 21) %>%
      bind_cols(s_america$xy) %>%
      mutate(region = "south america")
cpr_amazon <- cpr_rand_test(
      amazon$comm, amazon$tree, 
      null_model = "curveball", metrics = "pd", n_reps = 21) %>%
      bind_cols(amazon$xy) %>%
      mutate(region = "amazon")
cpr <- bind_rows(cpr_americas, cpr_s_america, cpr_amazon) %>%
      mutate(region = factor(region, levels = unique(region))) %>%
      mutate(pd_sig = case_when(pd_obs_p_upper > .95 ~ "high",
                                pd_obs_p_upper < .05 ~ "low",
                                TRUE ~ "no"))

# maps
cpr %>%
      carto("pd_sig") +
      facet_wrap(~ region, nrow = 1) +
      scale_fill_manual(values = c("red", "blue", "gray"))

cpr %>%
      filter(x %in% x[region == "amazon"],
             y %in% y[region == "amazon"]) %>%
      ggplot(aes(x, y, fill = pd_sig)) +
      geom_raster() +
      facet_wrap(~ region, nrow = 1) +
      theme_void() + theme(legend.position = "top") + coord_fixed() +
      scale_fill_manual(values = c("red", "blue", "gray"))

Facets of phylogenetic diversity

Phylogenetic branch lengths can represent different macroevolutionary variables, each with different implications for biogeography and conservation. We’ve indirectly met some different “phylogenetic diversity facets” already, through RPD (which is computed in canaper). Let’s quickly explore them in a bit more detail.

ds <- flora

# function to normalize branch lengths
scale_edges <- function(tree){
      tree$edge.length <- tree$edge.length / sum(tree$edge.length)
      return(tree)
}

# load trees
chronogram <- read.nexus("../../data/thornhill_2017/thornhill_2017_chronogram.nex") %>%
      scale_edges()
phylogram <- read.nexus("../../data/thornhill_2017/thornhill_2017_phylogram.nex") %>%
      scale_edges()

# make cladogram
cladogram <- phylogram
cladogram$edge.length <- rep(1, length(cladogram$edge.length))
cladogram <- scale_edges(cladogram)

# plot trees
par(mfrow = c(1, 3))
phyplot(chronogram, chronogram$edge.length,
        type = "fan", show.tip.label = F, main = "time")
phyplot(phylogram, phylogram$edge.length,
        type = "fan", show.tip.label = F, main = "divergence")
phyplot(cladogram, cladogram$edge.length,
        type = "fan", show.tip.label = F, main = "cladogenesis")

dev.off()
## null device 
##           1
# format as sparse matrix for phyloregion package
scomm <- ds$comm %>% Matrix(sparse = TRUE)

# PD for the three trees
div <- ds$xy %>%
      mutate(divergence = PD(scomm, phylogram),
             time = PD(scomm, chronogram),
             diversification = PD(scomm, cladogram))

# visualize
div %>%
      pivot_longer(c(divergence, time, diversification), 
                   "facet", "value") %>%
      group_by(facet) %>%
      mutate(rank = rank(value)) %>%
      carto("rank") + 
      facet_wrap(~facet, nrow = 1)
div %>%
      mutate(rpd = time / diversification) %>% 
      carto("rpd")
div %>%
      mutate(dynamism = divergence / time) %>%
      carto("dynamism")

Exercises:

  • Verify that rpd as returned by the cpr_rand_test() function using a chronogram is identical to the ratio of chronogram PD to cladogram PD calculated immediately above.
  • Make a map of the ratio of phylogram PD to cladogram PD, which measures the relative contribution of anagenesis versus cladogenesis to a community’s history. Where is this highest and lowest?
# RPD from cpr_rand_test() == chronogram PD / cladogram PD
cpr <- cpr_rand_test(ds$comm, ds$tree, null_model = "curveball", n_reps = 1)
pd_chron <- PD(scomm, chronogram)
pd_clad <- PD(scomm, cladogram)
plot(div$time/div$diversification, cpr$rpd_obs)

# phylogram PD / cladogram PD
div %>%
      mutate(genesis = divergence / diversification) %>%
      carto("genesis")

# phylogram PE / cladogram PE
ds$xy %>%
      mutate(divergence = phylo_endemism(scomm, phylogram),
             time = phylo_endemism(scomm, chronogram),
             diversification = phylo_endemism(scomm, cladogram)) %>%
      mutate(genesis = divergence / diversification) %>%
      carto("genesis")

BETA DIVERSITY

Analyses of beta diversity focus on compositional differences among local communities. Phylogenetic versions of these analyses represent differences between communities as the difference in the portion of the phylogenetic tree that’s found in each site.

Turnover

We’ll begin by creating a pairwise turnover matrix with a row and column for each site. We can compare this to species-based turnover and geographic distance. Further exploration of these turnover patterns could use a modeling approach like GDM (Ferrier 2007) to understand how factors like environmental differences explain phylogenetic turnover, but we will not cover that here.

ds <- mammals

# construct dissimlarity matrices
phy_beta <- ds$comm %>% # phylo sorensen's
      Matrix(sparse = TRUE) %>%
      phylobeta(ds$tree) %>%
      pluck("phylo.beta.sor") 
str(phy_beta)
##  'dist' Named num [1:3373503] 0.123 0.23 0.23 0.23 0.23 ...
##  - attr(*, "Labels")= chr [1:2598] "cell_1" "cell_2" "cell_3" "cell_4" ...
##  - attr(*, "Size")= int 2598
##  - attr(*, "Diag")= logi FALSE
##  - attr(*, "Upper")= logi FALSE
sp_beta <- ds$comm %>% # species sorensen's
      betapart::beta.pair() %>%
      pluck("beta.sor")
geo_dist <- ds$xy %>% # geographic distance
      geosphere::distm()

# plot geographic distance against phylogenetic turnover
ss <- sample(length(geo_dist), 10000)
tibble(distance = geo_dist[ss],
       phylo_sor = as.matrix(phy_beta)[ss],
       species_sor = as.matrix(sp_beta)[ss]) %>%
      pivot_longer(-distance, "stat", "value") %>%
      ggplot(aes(distance, value, color = stat)) +
      geom_point(size = .5) +
      geom_smooth(aes(group = stat), color = "black")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Regionalization

We can also use turnover data to make maps of phylogenetic similarity. Let’s look at two approaches: a discrete classification of phylogenetic regions using a cluster analysis from the phyloregions package, and a continuous ordination map that plots similar sites in similar colors:

# phylogenetic regionalization (spatial cluster analysis)
k <- 25
phyloregion(phy_beta, k = k) %>%
      pluck("membership") %>%
      bind_cols(ds$xy) %>%
      mutate(cluster = factor(cluster)) %>%
      carto("cluster") +
      scale_fill_manual(values = unname(pals::alphabet(k)))

# nmds of regions
    reg <- phyloregion(phy_beta, k = k)
    metaMDS(reg$region.dist, k = 2)$points %>%
      as.data.frame() %>%
      ggplot(aes(MDS1, MDS2, color = factor(1:k), label = 1:k)) +
      geom_label(fontface = "bold") +
      scale_color_manual(values = unname(pals::alphabet(k))) +
      theme(legend.position = "none")
## Run 0 stress 0.2018687 
## Run 1 stress 0.1991909 
## ... New best solution
## ... Procrustes: rmse 0.1412463  max resid 0.3582649 
## Run 2 stress 0.2018927 
## Run 3 stress 0.2017088 
## Run 4 stress 0.213275 
## Run 5 stress 0.208749 
## Run 6 stress 0.2173121 
## Run 7 stress 0.1893562 
## ... New best solution
## ... Procrustes: rmse 0.08448508  max resid 0.2859639 
## Run 8 stress 0.1975499 
## Run 9 stress 0.199811 
## Run 10 stress 0.2302428 
## Run 11 stress 0.2018688 
## Run 12 stress 0.224134 
## Run 13 stress 0.2169383 
## Run 14 stress 0.2168547 
## Run 15 stress 0.2082413 
## Run 16 stress 0.1997677 
## Run 17 stress 0.1874266 
## ... New best solution
## ... Procrustes: rmse 0.04206265  max resid 0.1656177 
## Run 18 stress 0.1972949 
## Run 19 stress 0.2061148 
## Run 20 stress 0.1975499 
## *** No convergence -- monoMDS stopping criteria:
##     17: stress ratio > sratmax
##      3: scale factor of the gradient < sfgrmin

    # unpacking and visualizing regional hclust
    library(dendextend)
## Registered S3 method overwritten by 'dendextend':
##   method     from 
##   rev.hclust vegan
## 
## ---------------------
## Welcome to dendextend version 1.16.0
## Type citation('dendextend') for how to cite the package.
## 
## Type browseVignettes(package = 'dendextend') for the package vignette.
## The github page is: https://github.com/talgalili/dendextend/
## 
## Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
## You may ask questions at stackoverflow, use the r and dendextend tags: 
##   https://stackoverflow.com/questions/tagged/dendextend
## 
##  To suppress this message use:  suppressPackageStartupMessages(library(dendextend))
## ---------------------
## 
## Attaching package: 'dendextend'
## The following object is masked from 'package:permute':
## 
##     shuffle
## The following object is masked from 'package:raster':
## 
##     rotate
## The following object is masked from 'package:terra':
## 
##     rotate
## The following objects are masked from 'package:ape':
## 
##     ladderize, rotate
## The following object is masked from 'package:stats':
## 
##     cutree
    hc <- phy_beta %>% hclust(method = "average")
    clust <- phyloregion(phy_beta, k = k) %>%
          pluck("membership")
    clrs <- unname(pals::alphabet(k))[clust$cluster]
    dend <- as.dendrogram(hc)
    o <- order.dendrogram(dend)
    dend <- assign_values_to_leaves_edgePar(
          dend, value = clrs[o], edgePar = "col")
    plot(dend, leaflab = "none")

    circlize_dendrogram(dend, labels = F, 
                        dend_track_height = .9, 
                        labels_track_height = 0)
## Loading required namespace: circlize

# RGB visualization
# ord <- vegan::metaMDS(phy_beta, k = 3, trymax = 50)$points
# saveRDS(ord$points, "../../data/cache/mammals_ord.rds")
ord <- readRDS("../../data/cache/mammals_ord.rds")
ord %>%
      as.data.frame() %>%
      mutate_all(function(x) rank(x)/max(rank(x))) %>%
      bind_cols(ds$xy) %>%
      mutate(color = rgb(MDS1, MDS2, MDS3)) %>%
      carto("color") +
      scale_fill_identity()

Exercises:

  • The above demonstration used Sorensen’s distance as a measure of phylogenetic community difference. Try things with Jaccard’s distance. How do the results differ?
  • Experiment with different numbers of clusters and different cluster methods in the phyloregion() function, to explore how robust regional boundaries are to these parameter choices.

CONSERVATION PRIORITIZATION

Phylogenetic diversity is useful for basic research, but from the beginning it has also been proposed as an applied conservation tool. Let’s use it to run a conservation prioritization analysis. We’ll focus on the California vascular flora dataset for this exercise, which has been used for conservation prioritization in the past (Kling et al. 2018).

Optimal reserve design is a major computational challenge, due to the extremely large number of sets of sites that could comprise a protected area network. Our computational workhorse will be the prioritizr library, a very flexible and powerful conservation planning toolkit that uses linear solvers to identify optimal conservation solutions.

First we need to get our data formatted as raster objects. We’ll also load an additional data layer of current protected areas, derived from Kling et al. (2018).

ds <- flora

# conservation features
features <- ds$xy %>%
      bind_cols(ds$comm) %>%
      rasterFromXYZ()

# spatial planning units
units <- features[[1]] %>% 
      reclassify(c(-Inf, Inf, 1))

# current protected areas
protected <- raster("../../data/thornhill_2017/kling_2018_protection_status.tif")
# plot(protected)
protected <- protected > .66
prot_area <- sum(protected[], na.rm = T)

# harmonize spatial coverage
overlap <- protected + units
protected <- protected %>% crop(overlap) %>% mask(overlap)
units <- units %>% crop(overlap) %>% mask(overlap)
features <- features %>% crop(overlap) %>% mask(overlap)
crs(protected) <- crs(units) <- crs(features)

The prioritizr library has some native functionality for incorporating phylogenies into conservation prioritization, so let’s start there.

# compare existing protected area network to 
# an optimal network of the same size
prob <- problem(units, features) %>%
      add_max_phylo_div_objective(budget = prot_area,
                                  tree = ds$tree) %>%
      add_relative_targets(0.1) %>%
      add_binary_decisions() %>%
      add_rsymphony_solver(gap = 0)
soln <- solve(prob)
stack(soln, protected) %>% 
      setNames(c("optimal", "existing")) %>%
      plot(col = c("gray", "orange"))

# identify the highest-priority sites to add to existing protected areas
prob <- prob %>%
      add_max_phylo_div_objective(budget = prot_area + 10,
                                  tree = ds$tree) %>%
      add_locked_in_constraints(locked_in = protected)
## Warning in res(x, ...): overwriting previously defined objective
soln <- solve(prob)
plot(soln + protected, 
     col = c("gray", "orange", "darkgreen"))

While it incorporates phylogeny, this method is still “tip-centric” in the sense that its objective is to maximize the phylogenetic diversity of the tips protected across the reserve network – the targets are still the tips, rather than taxa at all levels. A subtly different alternative that’s arguably more consistent with spatial phylogenetic metrics like PE is to build out a phylogenetic community matrix, and weight each taxon by its branch length during prioritization. We can see this yields a somewhat different result:

# build feature set representing all clades, not just tips
pcomm <- build_clade_ranges(ds$tree, ds$comm)
pfeatures <- ds$xy %>%
      bind_cols(pcomm) %>%
      rasterFromXYZ()
## New names:
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`
## • `` -> `...14`
## • `` -> `...15`
## • `` -> `...16`
## • `` -> `...17`
## • `` -> `...18`
## • `` -> `...19`
## • `` -> `...20`
## • `` -> `...21`
## • `` -> `...22`
## • `` -> `...23`
## • `` -> `...24`
## • `` -> `...25`
## • `` -> `...26`
## • `` -> `...27`
## • `` -> `...28`
## • `` -> `...29`
## • `` -> `...30`
## • `` -> `...31`
## • `` -> `...32`
## • `` -> `...33`
## • `` -> `...34`
## • `` -> `...35`
## • `` -> `...36`
## • `` -> `...37`
## • `` -> `...38`
## • `` -> `...39`
## • `` -> `...40`
## • `` -> `...41`
## • `` -> `...42`
## • `` -> `...43`
## • `` -> `...44`
## • `` -> `...45`
## • `` -> `...46`
## • `` -> `...47`
## • `` -> `...48`
## • `` -> `...49`
## • `` -> `...50`
## • `` -> `...51`
## • `` -> `...52`
## • `` -> `...53`
## • `` -> `...54`
## • `` -> `...55`
## • `` -> `...56`
## • `` -> `...57`
## • `` -> `...58`
## • `` -> `...59`
## • `` -> `...60`
## • `` -> `...61`
## • `` -> `...62`
## • `` -> `...63`
## • `` -> `...64`
## • `` -> `...65`
## • `` -> `...66`
## • `` -> `...67`
## • `` -> `...68`
## • `` -> `...69`
## • `` -> `...70`
## • `` -> `...71`
## • `` -> `...72`
## • `` -> `...73`
## • `` -> `...74`
## • `` -> `...75`
## • `` -> `...76`
## • `` -> `...77`
## • `` -> `...78`
## • `` -> `...79`
## • `` -> `...80`
## • `` -> `...81`
## • `` -> `...82`
## • `` -> `...83`
## • `` -> `...84`
## • `` -> `...85`
## • `` -> `...86`
## • `` -> `...87`
## • `` -> `...88`
## • `` -> `...89`
## • `` -> `...90`
## • `` -> `...91`
## • `` -> `...92`
## • `` -> `...93`
## • `` -> `...94`
## • `` -> `...95`
## • `` -> `...96`
## • `` -> `...97`
## • `` -> `...98`
## • `` -> `...99`
## • `` -> `...100`
## • `` -> `...101`
## • `` -> `...102`
## • `` -> `...103`
## • `` -> `...104`
## • `` -> `...105`
## • `` -> `...106`
## • `` -> `...107`
## • `` -> `...108`
## • `` -> `...109`
## • `` -> `...110`
## • `` -> `...111`
## • `` -> `...112`
## • `` -> `...113`
## • `` -> `...114`
## • `` -> `...115`
## • `` -> `...116`
## • `` -> `...117`
## • `` -> `...118`
## • `` -> `...119`
## • `` -> `...120`
## • `` -> `...121`
## • `` -> `...122`
## • `` -> `...123`
## • `` -> `...124`
## • `` -> `...125`
## • `` -> `...126`
## • `` -> `...127`
## • `` -> `...128`
## • `` -> `...129`
## • `` -> `...130`
## • `` -> `...131`
## • `` -> `...132`
## • `` -> `...133`
## • `` -> `...134`
## • `` -> `...135`
## • `` -> `...136`
## • `` -> `...137`
## • `` -> `...138`
## • `` -> `...139`
## • `` -> `...140`
## • `` -> `...141`
## • `` -> `...142`
## • `` -> `...143`
## • `` -> `...144`
## • `` -> `...145`
## • `` -> `...146`
## • `` -> `...147`
## • `` -> `...148`
## • `` -> `...149`
## • `` -> `...150`
## • `` -> `...151`
## • `` -> `...152`
## • `` -> `...153`
## • `` -> `...154`
## • `` -> `...155`
## • `` -> `...156`
## • `` -> `...157`
## • `` -> `...158`
## • `` -> `...159`
## • `` -> `...160`
## • `` -> `...161`
## • `` -> `...162`
## • `` -> `...163`
## • `` -> `...164`
## • `` -> `...165`
## • `` -> `...166`
## • `` -> `...167`
## • `` -> `...168`
## • `` -> `...169`
## • `` -> `...170`
## • `` -> `...171`
## • `` -> `...172`
## • `` -> `...173`
## • `` -> `...174`
## • `` -> `...175`
## • `` -> `...176`
## • `` -> `...177`
## • `` -> `...178`
## • `` -> `...179`
## • `` -> `...180`
## • `` -> `...181`
## • `` -> `...182`
## • `` -> `...183`
## • `` -> `...184`
## • `` -> `...185`
## • `` -> `...186`
## • `` -> `...187`
## • `` -> `...188`
## • `` -> `...189`
## • `` -> `...190`
## • `` -> `...191`
## • `` -> `...192`
## • `` -> `...193`
## • `` -> `...194`
## • `` -> `...195`
## • `` -> `...196`
## • `` -> `...197`
## • `` -> `...198`
## • `` -> `...199`
## • `` -> `...200`
## • `` -> `...201`
## • `` -> `...202`
## • `` -> `...203`
## • `` -> `...204`
## • `` -> `...205`
## • `` -> `...206`
## • `` -> `...207`
## • `` -> `...208`
## • `` -> `...209`
## • `` -> `...210`
## • `` -> `...211`
## • `` -> `...212`
## • `` -> `...213`
## • `` -> `...214`
## • `` -> `...215`
## • `` -> `...216`
## • `` -> `...217`
## • `` -> `...218`
## • `` -> `...219`
## • `` -> `...220`
## • `` -> `...221`
## • `` -> `...222`
## • `` -> `...223`
## • `` -> `...224`
## • `` -> `...225`
## • `` -> `...226`
## • `` -> `...227`
## • `` -> `...228`
## • `` -> `...229`
## • `` -> `...230`
## • `` -> `...231`
## • `` -> `...232`
## • `` -> `...233`
## • `` -> `...234`
## • `` -> `...235`
## • `` -> `...236`
## • `` -> `...237`
## • `` -> `...238`
## • `` -> `...239`
## • `` -> `...240`
## • `` -> `...241`
## • `` -> `...242`
## • `` -> `...243`
## • `` -> `...244`
## • `` -> `...245`
## • `` -> `...246`
## • `` -> `...247`
## • `` -> `...248`
## • `` -> `...249`
## • `` -> `...250`
## • `` -> `...251`
## • `` -> `...252`
## • `` -> `...253`
## • `` -> `...254`
## • `` -> `...255`
## • `` -> `...256`
## • `` -> `...257`
## • `` -> `...258`
## • `` -> `...259`
## • `` -> `...260`
## • `` -> `...261`
## • `` -> `...262`
## • `` -> `...263`
## • `` -> `...264`
## • `` -> `...265`
## • `` -> `...266`
## • `` -> `...267`
## • `` -> `...268`
## • `` -> `...269`
## • `` -> `...270`
## • `` -> `...271`
## • `` -> `...272`
## • `` -> `...273`
## • `` -> `...274`
## • `` -> `...275`
## • `` -> `...276`
## • `` -> `...277`
## • `` -> `...278`
## • `` -> `...279`
## • `` -> `...280`
## • `` -> `...281`
## • `` -> `...282`
## • `` -> `...283`
## • `` -> `...284`
## • `` -> `...285`
## • `` -> `...286`
## • `` -> `...287`
## • `` -> `...288`
## • `` -> `...289`
## • `` -> `...290`
## • `` -> `...291`
## • `` -> `...292`
## • `` -> `...293`
## • `` -> `...294`
## • `` -> `...295`
## • `` -> `...296`
## • `` -> `...297`
## • `` -> `...298`
## • `` -> `...299`
## • `` -> `...300`
## • `` -> `...301`
## • `` -> `...302`
## • `` -> `...303`
## • `` -> `...304`
## • `` -> `...305`
## • `` -> `...306`
## • `` -> `...307`
## • `` -> `...308`
## • `` -> `...309`
## • `` -> `...310`
## • `` -> `...311`
## • `` -> `...312`
## • `` -> `...313`
## • `` -> `...314`
## • `` -> `...315`
## • `` -> `...316`
## • `` -> `...317`
## • `` -> `...318`
## • `` -> `...319`
## • `` -> `...320`
## • `` -> `...321`
## • `` -> `...322`
## • `` -> `...323`
## • `` -> `...324`
## • `` -> `...325`
## • `` -> `...326`
## • `` -> `...327`
## • `` -> `...328`
## • `` -> `...329`
## • `` -> `...330`
## • `` -> `...331`
## • `` -> `...332`
## • `` -> `...333`
## • `` -> `...334`
## • `` -> `...335`
## • `` -> `...336`
## • `` -> `...337`
## • `` -> `...338`
## • `` -> `...339`
## • `` -> `...340`
## • `` -> `...341`
## • `` -> `...342`
## • `` -> `...343`
## • `` -> `...344`
## • `` -> `...345`
## • `` -> `...346`
## • `` -> `...347`
## • `` -> `...348`
## • `` -> `...349`
## • `` -> `...350`
## • `` -> `...351`
## • `` -> `...352`
## • `` -> `...353`
## • `` -> `...354`
## • `` -> `...355`
## • `` -> `...356`
## • `` -> `...357`
## • `` -> `...358`
## • `` -> `...359`
## • `` -> `...360`
## • `` -> `...361`
## • `` -> `...362`
## • `` -> `...363`
## • `` -> `...364`
## • `` -> `...365`
## • `` -> `...366`
## • `` -> `...367`
## • `` -> `...368`
## • `` -> `...369`
## • `` -> `...370`
## • `` -> `...371`
## • `` -> `...372`
## • `` -> `...373`
## • `` -> `...374`
## • `` -> `...375`
## • `` -> `...376`
## • `` -> `...377`
## • `` -> `...378`
## • `` -> `...379`
## • `` -> `...380`
## • `` -> `...381`
## • `` -> `...382`
## • `` -> `...383`
## • `` -> `...384`
## • `` -> `...385`
## • `` -> `...386`
## • `` -> `...387`
## • `` -> `...388`
## • `` -> `...389`
## • `` -> `...390`
## • `` -> `...391`
## • `` -> `...392`
## • `` -> `...393`
## • `` -> `...394`
## • `` -> `...395`
## • `` -> `...396`
## • `` -> `...397`
## • `` -> `...398`
## • `` -> `...399`
## • `` -> `...400`
## • `` -> `...401`
## • `` -> `...402`
## • `` -> `...403`
## • `` -> `...404`
## • `` -> `...405`
## • `` -> `...406`
## • `` -> `...407`
## • `` -> `...408`
## • `` -> `...409`
## • `` -> `...410`
## • `` -> `...411`
## • `` -> `...412`
## • `` -> `...413`
## • `` -> `...414`
## • `` -> `...415`
## • `` -> `...416`
## • `` -> `...417`
## • `` -> `...418`
## • `` -> `...419`
## • `` -> `...420`
## • `` -> `...421`
## • `` -> `...422`
## • `` -> `...423`
## • `` -> `...424`
## • `` -> `...425`
## • `` -> `...426`
## • `` -> `...427`
## • `` -> `...428`
## • `` -> `...429`
## • `` -> `...430`
## • `` -> `...431`
## • `` -> `...432`
## • `` -> `...433`
## • `` -> `...434`
## • `` -> `...435`
## • `` -> `...436`
## • `` -> `...437`
## • `` -> `...438`
## • `` -> `...439`
## • `` -> `...440`
## • `` -> `...441`
## • `` -> `...442`
## • `` -> `...443`
## • `` -> `...444`
## • `` -> `...445`
## • `` -> `...446`
## • `` -> `...447`
## • `` -> `...448`
## • `` -> `...449`
## • `` -> `...450`
## • `` -> `...451`
## • `` -> `...452`
## • `` -> `...453`
## • `` -> `...454`
## • `` -> `...455`
## • `` -> `...456`
## • `` -> `...457`
## • `` -> `...458`
## • `` -> `...459`
## • `` -> `...460`
## • `` -> `...461`
## • `` -> `...462`
## • `` -> `...463`
## • `` -> `...464`
## • `` -> `...465`
## • `` -> `...466`
## • `` -> `...467`
## • `` -> `...468`
## • `` -> `...469`
## • `` -> `...470`
## • `` -> `...471`
## • `` -> `...472`
## • `` -> `...473`
## • `` -> `...474`
## • `` -> `...475`
## • `` -> `...476`
## • `` -> `...477`
## • `` -> `...478`
## • `` -> `...479`
## • `` -> `...480`
## • `` -> `...481`
## • `` -> `...482`
## • `` -> `...483`
## • `` -> `...484`
## • `` -> `...485`
## • `` -> `...486`
## • `` -> `...487`
## • `` -> `...488`
## • `` -> `...489`
## • `` -> `...490`
## • `` -> `...491`
## • `` -> `...492`
## • `` -> `...493`
## • `` -> `...494`
## • `` -> `...495`
## • `` -> `...496`
## • `` -> `...497`
## • `` -> `...498`
## • `` -> `...499`
## • `` -> `...500`
## • `` -> `...501`
## • `` -> `...502`
## • `` -> `...503`
## • `` -> `...504`
## • `` -> `...505`
## • `` -> `...506`
## • `` -> `...507`
## • `` -> `...508`
## • `` -> `...509`
## • `` -> `...510`
## • `` -> `...511`
## • `` -> `...512`
## • `` -> `...513`
## • `` -> `...514`
## • `` -> `...515`
## • `` -> `...516`
## • `` -> `...517`
## • `` -> `...518`
## • `` -> `...519`
## • `` -> `...520`
## • `` -> `...521`
## • `` -> `...522`
## • `` -> `...523`
## • `` -> `...524`
## • `` -> `...525`
## • `` -> `...526`
## • `` -> `...527`
## • `` -> `...528`
## • `` -> `...529`
## • `` -> `...530`
## • `` -> `...531`
## • `` -> `...532`
## • `` -> `...533`
## • `` -> `...534`
## • `` -> `...535`
## • `` -> `...536`
## • `` -> `...537`
## • `` -> `...538`
## • `` -> `...539`
## • `` -> `...540`
## • `` -> `...541`
## • `` -> `...542`
## • `` -> `...543`
## • `` -> `...544`
## • `` -> `...545`
## • `` -> `...546`
## • `` -> `...547`
## • `` -> `...548`
## • `` -> `...549`
## • `` -> `...550`
## • `` -> `...551`
## • `` -> `...552`
## • `` -> `...553`
## • `` -> `...554`
## • `` -> `...555`
## • `` -> `...556`
## • `` -> `...557`
## • `` -> `...558`
## • `` -> `...559`
## • `` -> `...560`
## • `` -> `...561`
## • `` -> `...562`
## • `` -> `...563`
## • `` -> `...564`
## • `` -> `...565`
## • `` -> `...566`
## • `` -> `...567`
## • `` -> `...568`
## • `` -> `...569`
## • `` -> `...570`
## • `` -> `...571`
## • `` -> `...572`
## • `` -> `...573`
## • `` -> `...574`
## • `` -> `...575`
## • `` -> `...576`
## • `` -> `...577`
## • `` -> `...578`
## • `` -> `...579`
## • `` -> `...580`
## • `` -> `...581`
## • `` -> `...582`
## • `` -> `...583`
## • `` -> `...584`
## • `` -> `...585`
## • `` -> `...586`
## • `` -> `...587`
## • `` -> `...588`
## • `` -> `...589`
## • `` -> `...590`
## • `` -> `...591`
## • `` -> `...592`
## • `` -> `...593`
## • `` -> `...594`
## • `` -> `...595`
## • `` -> `...596`
## • `` -> `...597`
## • `` -> `...598`
## • `` -> `...599`
## • `` -> `...600`
## • `` -> `...601`
## • `` -> `...602`
## • `` -> `...603`
## • `` -> `...604`
## • `` -> `...605`
## • `` -> `...606`
## • `` -> `...607`
## • `` -> `...608`
## • `` -> `...609`
## • `` -> `...610`
## • `` -> `...611`
## • `` -> `...612`
## • `` -> `...613`
## • `` -> `...614`
## • `` -> `...615`
## • `` -> `...616`
## • `` -> `...617`
## • `` -> `...618`
## • `` -> `...619`
## • `` -> `...620`
## • `` -> `...621`
## • `` -> `...622`
## • `` -> `...623`
## • `` -> `...624`
## • `` -> `...625`
## • `` -> `...626`
## • `` -> `...627`
## • `` -> `...628`
## • `` -> `...629`
## • `` -> `...630`
## • `` -> `...631`
## • `` -> `...632`
## • `` -> `...633`
## • `` -> `...634`
## • `` -> `...635`
## • `` -> `...636`
## • `` -> `...637`
## • `` -> `...638`
## • `` -> `...639`
## • `` -> `...640`
## • `` -> `...641`
## • `` -> `...642`
## • `` -> `...643`
## • `` -> `...644`
## • `` -> `...645`
## • `` -> `...646`
## • `` -> `...647`
## • `` -> `...648`
## • `` -> `...649`
## • `` -> `...650`
## • `` -> `...651`
## • `` -> `...652`
## • `` -> `...653`
## • `` -> `...654`
## • `` -> `...655`
## • `` -> `...656`
## • `` -> `...657`
## • `` -> `...658`
## • `` -> `...659`
## • `` -> `...660`
## • `` -> `...661`
## • `` -> `...662`
## • `` -> `...663`
## • `` -> `...664`
## • `` -> `...665`
## • `` -> `...666`
## • `` -> `...667`
## • `` -> `...668`
## • `` -> `...669`
## • `` -> `...670`
## • `` -> `...671`
## • `` -> `...672`
## • `` -> `...673`
## • `` -> `...674`
## • `` -> `...675`
## • `` -> `...676`
## • `` -> `...677`
## • `` -> `...678`
## • `` -> `...679`
## • `` -> `...680`
## • `` -> `...681`
## • `` -> `...682`
## • `` -> `...683`
## • `` -> `...684`
## • `` -> `...685`
## • `` -> `...686`
## • `` -> `...687`
## • `` -> `...688`
## • `` -> `...689`
## • `` -> `...690`
## • `` -> `...691`
## • `` -> `...692`
## • `` -> `...693`
## • `` -> `...694`
## • `` -> `...695`
## • `` -> `...696`
## • `` -> `...697`
## • `` -> `...698`
## • `` -> `...699`
## • `` -> `...700`
## • `` -> `...701`
## • `` -> `...702`
## • `` -> `...703`
## • `` -> `...704`
## • `` -> `...705`
## • `` -> `...706`
## • `` -> `...707`
## • `` -> `...708`
## • `` -> `...709`
## • `` -> `...710`
## • `` -> `...711`
## • `` -> `...712`
## • `` -> `...713`
## • `` -> `...714`
## • `` -> `...715`
## • `` -> `...716`
## • `` -> `...717`
## • `` -> `...718`
## • `` -> `...719`
## • `` -> `...720`
## • `` -> `...721`
## • `` -> `...722`
## • `` -> `...723`
## • `` -> `...724`
## • `` -> `...725`
## • `` -> `...726`
## • `` -> `...727`
## • `` -> `...728`
## • `` -> `...729`
## • `` -> `...730`
## • `` -> `...731`
## • `` -> `...732`
## • `` -> `...733`
## • `` -> `...734`
## • `` -> `...735`
## • `` -> `...736`
## • `` -> `...737`
## • `` -> `...738`
## • `` -> `...739`
## • `` -> `...740`
## • `` -> `...741`
## • `` -> `...742`
## • `` -> `...743`
## • `` -> `...744`
## • `` -> `...745`
## • `` -> `...746`
## • `` -> `...747`
## • `` -> `...748`
## • `` -> `...749`
## • `` -> `...750`
## • `` -> `...751`
## • `` -> `...752`
## • `` -> `...753`
## • `` -> `...754`
## • `` -> `...755`
## • `` -> `...756`
## • `` -> `...757`
## • `` -> `...758`
## • `` -> `...759`
## • `` -> `...760`
## • `` -> `...761`
## • `` -> `...762`
## • `` -> `...763`
## • `` -> `...764`
## • `` -> `...765`
## • `` -> `...766`
## • `` -> `...767`
## • `` -> `...768`
## • `` -> `...769`
## • `` -> `...770`
## • `` -> `...771`
## • `` -> `...772`
## • `` -> `...773`
## • `` -> `...774`
## • `` -> `...775`
## • `` -> `...776`
## • `` -> `...777`
## • `` -> `...778`
## • `` -> `...779`
## • `` -> `...780`
## • `` -> `...781`
## • `` -> `...782`
## • `` -> `...783`
## • `` -> `...784`
## • `` -> `...785`
## • `` -> `...786`
## • `` -> `...787`
## • `` -> `...788`
## • `` -> `...789`
## • `` -> `...790`
## • `` -> `...791`
## • `` -> `...792`
## • `` -> `...793`
## • `` -> `...794`
## • `` -> `...795`
## • `` -> `...796`
## • `` -> `...797`
## • `` -> `...798`
## • `` -> `...799`
## • `` -> `...800`
## • `` -> `...801`
## • `` -> `...802`
## • `` -> `...803`
## • `` -> `...804`
## • `` -> `...805`
## • `` -> `...806`
## • `` -> `...807`
## • `` -> `...808`
## • `` -> `...809`
## • `` -> `...810`
## • `` -> `...811`
## • `` -> `...812`
## • `` -> `...813`
## • `` -> `...814`
## • `` -> `...815`
## • `` -> `...816`
## • `` -> `...817`
## • `` -> `...818`
## • `` -> `...819`
## • `` -> `...820`
## • `` -> `...821`
## • `` -> `...822`
## • `` -> `...823`
## • `` -> `...824`
## • `` -> `...825`
## • `` -> `...826`
## • `` -> `...827`
## • `` -> `...828`
## • `` -> `...829`
## • `` -> `...830`
## • `` -> `...831`
## • `` -> `...832`
## • `` -> `...833`
## • `` -> `...834`
## • `` -> `...835`
## • `` -> `...836`
## • `` -> `...837`
## • `` -> `...838`
## • `` -> `...839`
## • `` -> `...840`
## • `` -> `...841`
## • `` -> `...842`
## • `` -> `...843`
## • `` -> `...844`
## • `` -> `...845`
## • `` -> `...846`
## • `` -> `...847`
## • `` -> `...848`
## • `` -> `...849`
## • `` -> `...850`
## • `` -> `...851`
## • `` -> `...852`
## • `` -> `...853`
## • `` -> `...854`
## • `` -> `...855`
## • `` -> `...856`
## • `` -> `...857`
## • `` -> `...858`
## • `` -> `...859`
## • `` -> `...860`
## • `` -> `...861`
## • `` -> `...862`
## • `` -> `...863`
## • `` -> `...864`
## • `` -> `...865`
## • `` -> `...866`
## • `` -> `...867`
## • `` -> `...868`
## • `` -> `...869`
## • `` -> `...870`
## • `` -> `...871`
## • `` -> `...872`
## • `` -> `...873`
## • `` -> `...874`
## • `` -> `...875`
## • `` -> `...876`
## • `` -> `...877`
## • `` -> `...878`
## • `` -> `...879`
## • `` -> `...880`
## • `` -> `...881`
## • `` -> `...882`
## • `` -> `...883`
## • `` -> `...884`
## • `` -> `...885`
## • `` -> `...886`
## • `` -> `...887`
## • `` -> `...888`
## • `` -> `...889`
## • `` -> `...890`
## • `` -> `...891`
## • `` -> `...892`
## • `` -> `...893`
## • `` -> `...894`
## • `` -> `...895`
## • `` -> `...896`
## • `` -> `...897`
## • `` -> `...898`
## • `` -> `...899`
## • `` -> `...900`
## • `` -> `...901`
## • `` -> `...902`
## • `` -> `...903`
## • `` -> `...904`
## • `` -> `...905`
## • `` -> `...906`
## • `` -> `...907`
## • `` -> `...908`
## • `` -> `...909`
## • `` -> `...910`
## • `` -> `...911`
## • `` -> `...912`
## • `` -> `...913`
## • `` -> `...914`
## • `` -> `...915`
## • `` -> `...916`
## • `` -> `...917`
## • `` -> `...918`
## • `` -> `...919`
## • `` -> `...920`
## • `` -> `...921`
## • `` -> `...922`
## • `` -> `...923`
## • `` -> `...924`
## • `` -> `...925`
## • `` -> `...926`
## • `` -> `...927`
## • `` -> `...928`
## • `` -> `...929`
## • `` -> `...930`
## • `` -> `...931`
## • `` -> `...932`
## • `` -> `...933`
## • `` -> `...934`
## • `` -> `...935`
## • `` -> `...936`
## • `` -> `...937`
## • `` -> `...938`
## • `` -> `...939`
## • `` -> `...940`
## • `` -> `...941`
## • `` -> `...942`
## • `` -> `...943`
## • `` -> `...944`
## • `` -> `...945`
## • `` -> `...946`
## • `` -> `...947`
## • `` -> `...948`
## • `` -> `...949`
## • `` -> `...950`
## • `` -> `...951`
## • `` -> `...952`
## • `` -> `...953`
## • `` -> `...954`
## • `` -> `...955`
## • `` -> `...956`
## • `` -> `...957`
## • `` -> `...958`
## • `` -> `...959`
## • `` -> `...960`
## • `` -> `...961`
## • `` -> `...962`
## • `` -> `...963`
## • `` -> `...964`
## • `` -> `...965`
## • `` -> `...966`
## • `` -> `...967`
## • `` -> `...968`
## • `` -> `...969`
## • `` -> `...970`
## • `` -> `...971`
## • `` -> `...972`
## • `` -> `...973`
## • `` -> `...974`
## • `` -> `...975`
## • `` -> `...976`
## • `` -> `...977`
## • `` -> `...978`
## • `` -> `...979`
## • `` -> `...980`
## • `` -> `...981`
## • `` -> `...982`
## • `` -> `...983`
## • `` -> `...984`
## • `` -> `...985`
## • `` -> `...986`
## • `` -> `...987`
## • `` -> `...988`
## • `` -> `...989`
## • `` -> `...990`
## • `` -> `...991`
## • `` -> `...992`
## • `` -> `...993`
## • `` -> `...994`
## • `` -> `...995`
## • `` -> `...996`
## • `` -> `...997`
## • `` -> `...998`
## • `` -> `...999`
## • `` -> `...1000`
## • `` -> `...1001`
## • `` -> `...1002`
## • `` -> `...1003`
## • `` -> `...1004`
## • `` -> `...1005`
## • `` -> `...1006`
## • `` -> `...1007`
## • `` -> `...1008`
## • `` -> `...1009`
## • `` -> `...1010`
## • `` -> `...1011`
## • `` -> `...1012`
## • `` -> `...1013`
## • `` -> `...1014`
## • `` -> `...1015`
## • `` -> `...1016`
## • `` -> `...1017`
## • `` -> `...1018`
## • `` -> `...1019`
## • `` -> `...1020`
## • `` -> `...1021`
## • `` -> `...1022`
## • `` -> `...1023`
## • `` -> `...1024`
## • `` -> `...1025`
## • `` -> `...1026`
## • `` -> `...1027`
## • `` -> `...1028`
## • `` -> `...1029`
## • `` -> `...1030`
## • `` -> `...1031`
## • `` -> `...1032`
## • `` -> `...1033`
## • `` -> `...1034`
## • `` -> `...1035`
## • `` -> `...1036`
## • `` -> `...1037`
## • `` -> `...1038`
## • `` -> `...1039`
## • `` -> `...1040`
## • `` -> `...1041`
## • `` -> `...1042`
## • `` -> `...1043`
## • `` -> `...1044`
## • `` -> `...1045`
## • `` -> `...1046`
## • `` -> `...1047`
## • `` -> `...1048`
## • `` -> `...1049`
## • `` -> `...1050`
## • `` -> `...1051`
## • `` -> `...1052`
## • `` -> `...1053`
## • `` -> `...1054`
## • `` -> `...1055`
## • `` -> `...1056`
## • `` -> `...1057`
## • `` -> `...1058`
## • `` -> `...1059`
## • `` -> `...1060`
## • `` -> `...1061`
## • `` -> `...1062`
## • `` -> `...1063`
## • `` -> `...1064`
## • `` -> `...1065`
## • `` -> `...1066`
## • `` -> `...1067`
## • `` -> `...1068`
## • `` -> `...1069`
## • `` -> `...1070`
## • `` -> `...1071`
## • `` -> `...1072`
## • `` -> `...1073`
## • `` -> `...1074`
## • `` -> `...1075`
## • `` -> `...1076`
## • `` -> `...1077`
## • `` -> `...1078`
## • `` -> `...1079`
## • `` -> `...1080`
## • `` -> `...1081`
## • `` -> `...1082`
## • `` -> `...1083`
## • `` -> `...1084`
## • `` -> `...1085`
## • `` -> `...1086`
## • `` -> `...1087`
## • `` -> `...1088`
## • `` -> `...1089`
## • `` -> `...1090`
## • `` -> `...1091`
## • `` -> `...1092`
## • `` -> `...1093`
## • `` -> `...1094`
## • `` -> `...1095`
## • `` -> `...1096`
## • `` -> `...1097`
## • `` -> `...1098`
## • `` -> `...1099`
## • `` -> `...1100`
## • `` -> `...1101`
## • `` -> `...1102`
## • `` -> `...1103`
## • `` -> `...1104`
## • `` -> `...1105`
## • `` -> `...1106`
## • `` -> `...1107`
## • `` -> `...1108`
## • `` -> `...1109`
## • `` -> `...1110`
## • `` -> `...1111`
## • `` -> `...1112`
## • `` -> `...1113`
## • `` -> `...1114`
## • `` -> `...1115`
## • `` -> `...1116`
## • `` -> `...1117`
## • `` -> `...1118`
## • `` -> `...1119`
## • `` -> `...1120`
## • `` -> `...1121`
## • `` -> `...1122`
## • `` -> `...1123`
## • `` -> `...1124`
## • `` -> `...1125`
## • `` -> `...1126`
## • `` -> `...1127`
## • `` -> `...1128`
## • `` -> `...1129`
## • `` -> `...1130`
## • `` -> `...1131`
## • `` -> `...1132`
## • `` -> `...1133`
## • `` -> `...1134`
## • `` -> `...1135`
## • `` -> `...1136`
## • `` -> `...1137`
## • `` -> `...1138`
## • `` -> `...1139`
## • `` -> `...1140`
## • `` -> `...1141`
## • `` -> `...1142`
## • `` -> `...1143`
## • `` -> `...1144`
## • `` -> `...1145`
## • `` -> `...1146`
## • `` -> `...1147`
## • `` -> `...1148`
## • `` -> `...1149`
## • `` -> `...1150`
## • `` -> `...1151`
## • `` -> `...1152`
## • `` -> `...1153`
## • `` -> `...1154`
## • `` -> `...1155`
## • `` -> `...1156`
## • `` -> `...1157`
## • `` -> `...1158`
## • `` -> `...1159`
## • `` -> `...1160`
## • `` -> `...1161`
## • `` -> `...1162`
## • `` -> `...1163`
## • `` -> `...1164`
## • `` -> `...1165`
## • `` -> `...1166`
## • `` -> `...1167`
## • `` -> `...1168`
## • `` -> `...1169`
## • `` -> `...1170`
## • `` -> `...1171`
## • `` -> `...1172`
## • `` -> `...1173`
## • `` -> `...1174`
## • `` -> `...1175`
## • `` -> `...1176`
## • `` -> `...1177`
## • `` -> `...1178`
## • `` -> `...1179`
## • `` -> `...1180`
## • `` -> `...1181`
## • `` -> `...1182`
## • `` -> `...1183`
## • `` -> `...1184`
## • `` -> `...1185`
## • `` -> `...1186`
## • `` -> `...1187`
## • `` -> `...1188`
## • `` -> `...1189`
## • `` -> `...1190`
## • `` -> `...1191`
## • `` -> `...1192`
## • `` -> `...1193`
## • `` -> `...1194`
## • `` -> `...1195`
## • `` -> `...1196`
## • `` -> `...1197`
## • `` -> `...1198`
## • `` -> `...1199`
## • `` -> `...1200`
## • `` -> `...1201`
## • `` -> `...1202`
## • `` -> `...1203`
## • `` -> `...1204`
## • `` -> `...1205`
## • `` -> `...1206`
## • `` -> `...1207`
## • `` -> `...1208`
## • `` -> `...1209`
## • `` -> `...1210`
## • `` -> `...1211`
## • `` -> `...1212`
## • `` -> `...1213`
## • `` -> `...1214`
## • `` -> `...1215`
## • `` -> `...1216`
## • `` -> `...1217`
## • `` -> `...1218`
## • `` -> `...1219`
## • `` -> `...1220`
## • `` -> `...1221`
## • `` -> `...1222`
## • `` -> `...1223`
## • `` -> `...1224`
## • `` -> `...1225`
## • `` -> `...1226`
## • `` -> `...1227`
## • `` -> `...1228`
## • `` -> `...1229`
## • `` -> `...1230`
## • `` -> `...1231`
## • `` -> `...1232`
## • `` -> `...1233`
## • `` -> `...1234`
## • `` -> `...1235`
## • `` -> `...1236`
## • `` -> `...1237`
## • `` -> `...1238`
## • `` -> `...1239`
## • `` -> `...1240`
## • `` -> `...1241`
## • `` -> `...1242`
## • `` -> `...1243`
## • `` -> `...1244`
## • `` -> `...1245`
## • `` -> `...1246`
## • `` -> `...1247`
## • `` -> `...1248`
## • `` -> `...1249`
## • `` -> `...1250`
## • `` -> `...1251`
## • `` -> `...1252`
## • `` -> `...1253`
## • `` -> `...1254`
## • `` -> `...1255`
## • `` -> `...1256`
## • `` -> `...1257`
## • `` -> `...1258`
## • `` -> `...1259`
## • `` -> `...1260`
## • `` -> `...1261`
## • `` -> `...1262`
## • `` -> `...1263`
## • `` -> `...1264`
## • `` -> `...1265`
## • `` -> `...1266`
## • `` -> `...1267`
## • `` -> `...1268`
## • `` -> `...1269`
## • `` -> `...1270`
## • `` -> `...1271`
## • `` -> `...1272`
## • `` -> `...1273`
## • `` -> `...1274`
## • `` -> `...1275`
## • `` -> `...1276`
## • `` -> `...1277`
## • `` -> `...1278`
## • `` -> `...1279`
## • `` -> `...1280`
## • `` -> `...1281`
## • `` -> `...1282`
## • `` -> `...1283`
## • `` -> `...1284`
## • `` -> `...1285`
## • `` -> `...1286`
## • `` -> `...1287`
## • `` -> `...1288`
## • `` -> `...1289`
## • `` -> `...1290`
## • `` -> `...1291`
## • `` -> `...1292`
## • `` -> `...1293`
## • `` -> `...1294`
## • `` -> `...1295`
## • `` -> `...1296`
## • `` -> `...1297`
## • `` -> `...1298`
## • `` -> `...1299`
## • `` -> `...1300`
## • `` -> `...1301`
## • `` -> `...1302`
## • `` -> `...1303`
## • `` -> `...1304`
## • `` -> `...1305`
## • `` -> `...1306`
## • `` -> `...1307`
## • `` -> `...1308`
## • `` -> `...1309`
## • `` -> `...1310`
## • `` -> `...1311`
## • `` -> `...1312`
## • `` -> `...1313`
## • `` -> `...1314`
## • `` -> `...1315`
## • `` -> `...1316`
## • `` -> `...1317`
## • `` -> `...1318`
## • `` -> `...1319`
## • `` -> `...1320`
## • `` -> `...1321`
## • `` -> `...1322`
## • `` -> `...1323`
## • `` -> `...1324`
## • `` -> `...1325`
## • `` -> `...1326`
## • `` -> `...1327`
## • `` -> `...1328`
## • `` -> `...1329`
## • `` -> `...1330`
## • `` -> `...1331`
## • `` -> `...1332`
## • `` -> `...1333`
## • `` -> `...1334`
## • `` -> `...1335`
## • `` -> `...1336`
## • `` -> `...1337`
## • `` -> `...1338`
## • `` -> `...1339`
## • `` -> `...1340`
## • `` -> `...1341`
## • `` -> `...1342`
## • `` -> `...1343`
## • `` -> `...1344`
## • `` -> `...1345`
## • `` -> `...1346`
## • `` -> `...1347`
## • `` -> `...1348`
## • `` -> `...1349`
## • `` -> `...1350`
## • `` -> `...1351`
## • `` -> `...1352`
## • `` -> `...1353`
## • `` -> `...1354`
## • `` -> `...1355`
## • `` -> `...1356`
## • `` -> `...1357`
## • `` -> `...1358`
## • `` -> `...1359`
## • `` -> `...1360`
## • `` -> `...1361`
## • `` -> `...1362`
## • `` -> `...1363`
## • `` -> `...1364`
## • `` -> `...1365`
## • `` -> `...1366`
## • `` -> `...1367`
## • `` -> `...1368`
## • `` -> `...1369`
## • `` -> `...1370`
## • `` -> `...1371`
## • `` -> `...1372`
## • `` -> `...1373`
## • `` -> `...1374`
## • `` -> `...1375`
## • `` -> `...1376`
## • `` -> `...1377`
## • `` -> `...1378`
## • `` -> `...1379`
## • `` -> `...1380`
## • `` -> `...1381`
## • `` -> `...1382`
## • `` -> `...1383`
## • `` -> `...1384`
## • `` -> `...1385`
## • `` -> `...1386`
## • `` -> `...1387`
## • `` -> `...1388`
## • `` -> `...1389`
## • `` -> `...1390`
## • `` -> `...1391`
## • `` -> `...1392`
## • `` -> `...1393`
## • `` -> `...1394`
## • `` -> `...1395`
## • `` -> `...1396`
## • `` -> `...1397`
## • `` -> `...1398`
## • `` -> `...1399`
## • `` -> `...1400`
## • `` -> `...1401`
## • `` -> `...1402`
## • `` -> `...1403`
## • `` -> `...1404`
## • `` -> `...1405`
## • `` -> `...1406`
## • `` -> `...1407`
## • `` -> `...1408`
## • `` -> `...1409`
## • `` -> `...1410`
## • `` -> `...1411`
## • `` -> `...1412`
## • `` -> `...1413`
## • `` -> `...1414`
## • `` -> `...1415`
## • `` -> `...1416`
## • `` -> `...1417`
## • `` -> `...1418`
## • `` -> `...1419`
## • `` -> `...1420`
## • `` -> `...1421`
## • `` -> `...1422`
## • `` -> `...1423`
## • `` -> `...1424`
## • `` -> `...1425`
## • `` -> `...1426`
## • `` -> `...1427`
## • `` -> `...1428`
## • `` -> `...1429`
## • `` -> `...1430`
## • `` -> `...1431`
## • `` -> `...1432`
## • `` -> `...1433`
## • `` -> `...1434`
## • `` -> `...1435`
## • `` -> `...1436`
## • `` -> `...1437`
## • `` -> `...1438`
## • `` -> `...1439`
## • `` -> `...1440`
## • `` -> `...1441`
## • `` -> `...1442`
## • `` -> `...1443`
## • `` -> `...1444`
## • `` -> `...1445`
## • `` -> `...1446`
## • `` -> `...1447`
## • `` -> `...1448`
## • `` -> `...1449`
## • `` -> `...1450`
## • `` -> `...1451`
## • `` -> `...1452`
## • `` -> `...1453`
## • `` -> `...1454`
## • `` -> `...1455`
## • `` -> `...1456`
## • `` -> `...1457`
## • `` -> `...1458`
## • `` -> `...1459`
## • `` -> `...1460`
## • `` -> `...1461`
## • `` -> `...1462`
## • `` -> `...1463`
## • `` -> `...1464`
## • `` -> `...1465`
## • `` -> `...1466`
## • `` -> `...1467`
## • `` -> `...1468`
## • `` -> `...1469`
## • `` -> `...1470`
## • `` -> `...1471`
## • `` -> `...1472`
## • `` -> `...1473`
## • `` -> `...1474`
## • `` -> `...1475`
## • `` -> `...1476`
## • `` -> `...1477`
## • `` -> `...1478`
## • `` -> `...1479`
## • `` -> `...1480`
## • `` -> `...1481`
## • `` -> `...1482`
## • `` -> `...1483`
## • `` -> `...1484`
## • `` -> `...1485`
## • `` -> `...1486`
## • `` -> `...1487`
## • `` -> `...1488`
## • `` -> `...1489`
## • `` -> `...1490`
## • `` -> `...1491`
## • `` -> `...1492`
## • `` -> `...1493`
## • `` -> `...1494`
## • `` -> `...1495`
## • `` -> `...1496`
## • `` -> `...1497`
## • `` -> `...1498`
## • `` -> `...1499`
## • `` -> `...1500`
## • `` -> `...1501`
## • `` -> `...1502`
## • `` -> `...1503`
## • `` -> `...1504`
## • `` -> `...1505`
## • `` -> `...1506`
## • `` -> `...1507`
## • `` -> `...1508`
## • `` -> `...1509`
## • `` -> `...1510`
## • `` -> `...1511`
## • `` -> `...1512`
## • `` -> `...1513`
## • `` -> `...1514`
## • `` -> `...1515`
## • `` -> `...1516`
## • `` -> `...1517`
## • `` -> `...1518`
## • `` -> `...1519`
## • `` -> `...1520`
## • `` -> `...1521`
## • `` -> `...1522`
## • `` -> `...1523`
## • `` -> `...1524`
## • `` -> `...1525`
## • `` -> `...1526`
## • `` -> `...1527`
## • `` -> `...1528`
## • `` -> `...1529`
## • `` -> `...1530`
## • `` -> `...1531`
## • `` -> `...1532`
## • `` -> `...1533`
## • `` -> `...1534`
## • `` -> `...1535`
## • `` -> `...1536`
## • `` -> `...1537`
## • `` -> `...1538`
## • `` -> `...1539`
## • `` -> `...1540`
## • `` -> `...1541`
## • `` -> `...1542`
## • `` -> `...1543`
## • `` -> `...1544`
## • `` -> `...1545`
## • `` -> `...1546`
## • `` -> `...1547`
## • `` -> `...1548`
## • `` -> `...1549`
## • `` -> `...1550`
## • `` -> `...1551`
## • `` -> `...1552`
## • `` -> `...1553`
## • `` -> `...1554`
## • `` -> `...1555`
## • `` -> `...1556`
## • `` -> `...1557`
## • `` -> `...1558`
## • `` -> `...1559`
## • `` -> `...1560`
## • `` -> `...1561`
## • `` -> `...1562`
## • `` -> `...1563`
## • `` -> `...1564`
## • `` -> `...1565`
## • `` -> `...1566`
## • `` -> `...1567`
## • `` -> `...1568`
## • `` -> `...1569`
## • `` -> `...1570`
## • `` -> `...1571`
## • `` -> `...1572`
## • `` -> `...1573`
## • `` -> `...1574`
## • `` -> `...1575`
## • `` -> `...1576`
## • `` -> `...1577`
## • `` -> `...1578`
## • `` -> `...1579`
## • `` -> `...1580`
## • `` -> `...1581`
## • `` -> `...1582`
## • `` -> `...1583`
## • `` -> `...1584`
## • `` -> `...1585`
## • `` -> `...1586`
## • `` -> `...1587`
## • `` -> `...1588`
## • `` -> `...1589`
## • `` -> `...1590`
## • `` -> `...1591`
## • `` -> `...1592`
## • `` -> `...1593`
## • `` -> `...1594`
## • `` -> `...1595`
## • `` -> `...1596`
## • `` -> `...1597`
## • `` -> `...1598`
## • `` -> `...1599`
## • `` -> `...1600`
## • `` -> `...1601`
## • `` -> `...1602`
## • `` -> `...1603`
## • `` -> `...1604`
## • `` -> `...1605`
## • `` -> `...1606`
## • `` -> `...1607`
## • `` -> `...1608`
## • `` -> `...1609`
## • `` -> `...1610`
## • `` -> `...1611`
## • `` -> `...1612`
## • `` -> `...1613`
## • `` -> `...1614`
## • `` -> `...1615`
## • `` -> `...1616`
## • `` -> `...1617`
## • `` -> `...1618`
## • `` -> `...1619`
## • `` -> `...1620`
## • `` -> `...1621`
## • `` -> `...1622`
## • `` -> `...1623`
## • `` -> `...1624`
## • `` -> `...1625`
## • `` -> `...1626`
## • `` -> `...1627`
## • `` -> `...1628`
## • `` -> `...1629`
## • `` -> `...1630`
## • `` -> `...1631`
## • `` -> `...1632`
## • `` -> `...1633`
## • `` -> `...1634`
## • `` -> `...1635`
## • `` -> `...1636`
## • `` -> `...1637`
## • `` -> `...1638`
## • `` -> `...1639`
## • `` -> `...1640`
## • `` -> `...1641`
## • `` -> `...1642`
## • `` -> `...1643`
## • `` -> `...1644`
## • `` -> `...1645`
## • `` -> `...1646`
## • `` -> `...1647`
## • `` -> `...1648`
## • `` -> `...1649`
## • `` -> `...1650`
## • `` -> `...1651`
## • `` -> `...1652`
## • `` -> `...1653`
## • `` -> `...1654`
## • `` -> `...1655`
## • `` -> `...1656`
## • `` -> `...1657`
## • `` -> `...1658`
## • `` -> `...1659`
## • `` -> `...1660`
## • `` -> `...1661`
## • `` -> `...1662`
## • `` -> `...1663`
## • `` -> `...1664`
## • `` -> `...1665`
## • `` -> `...1666`
## • `` -> `...1667`
## • `` -> `...1668`
## • `` -> `...1669`
## • `` -> `...1670`
## • `` -> `...1671`
## • `` -> `...1672`
## • `` -> `...1673`
## • `` -> `...1674`
## • `` -> `...1675`
## • `` -> `...1676`
## • `` -> `...1677`
## • `` -> `...1678`
## • `` -> `...1679`
## • `` -> `...1680`
## • `` -> `...1681`
## • `` -> `...1682`
## • `` -> `...1683`
## • `` -> `...1684`
## • `` -> `...1685`
## • `` -> `...1686`
## • `` -> `...1687`
## • `` -> `...1688`
## • `` -> `...1689`
## • `` -> `...1690`
## • `` -> `...1691`
## • `` -> `...1692`
## • `` -> `...1693`
## • `` -> `...1694`
## • `` -> `...1695`
## • `` -> `...1696`
## • `` -> `...1697`
## • `` -> `...1698`
## • `` -> `...1699`
## • `` -> `...1700`
## • `` -> `...1701`
## • `` -> `...1702`
## • `` -> `...1703`
## • `` -> `...1704`
## • `` -> `...1705`
## • `` -> `...1706`
## • `` -> `...1707`
## • `` -> `...1708`
## • `` -> `...1709`
## • `` -> `...1710`
## • `` -> `...1711`
## • `` -> `...1712`
## • `` -> `...1713`
## • `` -> `...1714`
## • `` -> `...1715`
## • `` -> `...1716`
## • `` -> `...1717`
## • `` -> `...1718`
## • `` -> `...1719`
## • `` -> `...1720`
## • `` -> `...1721`
## • `` -> `...1722`
## • `` -> `...1723`
## • `` -> `...1724`
## • `` -> `...1725`
## • `` -> `...1726`
## • `` -> `...1727`
## • `` -> `...1728`
## • `` -> `...1729`
## • `` -> `...1730`
## • `` -> `...1731`
## • `` -> `...1732`
## • `` -> `...1733`
## • `` -> `...1734`
## • `` -> `...1735`
## • `` -> `...1736`
## • `` -> `...1737`
## • `` -> `...1738`
## • `` -> `...1739`
## • `` -> `...1740`
## • `` -> `...1741`
## • `` -> `...1742`
## • `` -> `...1743`
## • `` -> `...1744`
## • `` -> `...1745`
## • `` -> `...1746`
## • `` -> `...1747`
## • `` -> `...1748`
## • `` -> `...1749`
## • `` -> `...1750`
## • `` -> `...1751`
## • `` -> `...1752`
## • `` -> `...1753`
## • `` -> `...1754`
## • `` -> `...1755`
## • `` -> `...1756`
## • `` -> `...1757`
## • `` -> `...1758`
## • `` -> `...1759`
## • `` -> `...1760`
## • `` -> `...1761`
## • `` -> `...1762`
## • `` -> `...1763`
## • `` -> `...1764`
## • `` -> `...1765`
## • `` -> `...1766`
## • `` -> `...1767`
## • `` -> `...1768`
## • `` -> `...1769`
## • `` -> `...1770`
## • `` -> `...1771`
## • `` -> `...1772`
## • `` -> `...1773`
## • `` -> `...1774`
## • `` -> `...1775`
## • `` -> `...1776`
## • `` -> `...1777`
## • `` -> `...1778`
## • `` -> `...1779`
## • `` -> `...1780`
## • `` -> `...1781`
## • `` -> `...1782`
## • `` -> `...1783`
## • `` -> `...1784`
## • `` -> `...1785`
## • `` -> `...1786`
## • `` -> `...1787`
## • `` -> `...1788`
## • `` -> `...1789`
## • `` -> `...1790`
## • `` -> `...1791`
## • `` -> `...1792`
## • `` -> `...1793`
## • `` -> `...1794`
## • `` -> `...1795`
## • `` -> `...1796`
## • `` -> `...1797`
## • `` -> `...1798`
## • `` -> `...1799`
## • `` -> `...1800`
## • `` -> `...1801`
## • `` -> `...1802`
## • `` -> `...1803`
## • `` -> `...1804`
## • `` -> `...1805`
## • `` -> `...1806`
## • `` -> `...1807`
## • `` -> `...1808`
## • `` -> `...1809`
## • `` -> `...1810`
## • `` -> `...1811`
## • `` -> `...1812`
## • `` -> `...1813`
## • `` -> `...1814`
## • `` -> `...1815`
## • `` -> `...1816`
## • `` -> `...1817`
## • `` -> `...1818`
## • `` -> `...1819`
## • `` -> `...1820`
## • `` -> `...1821`
## • `` -> `...1822`
## • `` -> `...1823`
## • `` -> `...1824`
## • `` -> `...1825`
## • `` -> `...1826`
## • `` -> `...1827`
## • `` -> `...1828`
## • `` -> `...1829`
## • `` -> `...1830`
## • `` -> `...1831`
## • `` -> `...1832`
## • `` -> `...1833`
## • `` -> `...1834`
## • `` -> `...1835`
## • `` -> `...1836`
## • `` -> `...1837`
## • `` -> `...1838`
## • `` -> `...1839`
## • `` -> `...1840`
## • `` -> `...1841`
## • `` -> `...1842`
## • `` -> `...1843`
## • `` -> `...1844`
## • `` -> `...1845`
## • `` -> `...1846`
## • `` -> `...1847`
## • `` -> `...1848`
## • `` -> `...1849`
## • `` -> `...1850`
## • `` -> `...1851`
## • `` -> `...1852`
## • `` -> `...1853`
## • `` -> `...1854`
## • `` -> `...1855`
## • `` -> `...1856`
## • `` -> `...1857`
## • `` -> `...1858`
## • `` -> `...1859`
## • `` -> `...1860`
## • `` -> `...1861`
## • `` -> `...1862`
## • `` -> `...1863`
## • `` -> `...1864`
## • `` -> `...1865`
## • `` -> `...1866`
## • `` -> `...1867`
## • `` -> `...1868`
## • `` -> `...1869`
## • `` -> `...1870`
## • `` -> `...1871`
## • `` -> `...1872`
## • `` -> `...1873`
## • `` -> `...1874`
## • `` -> `...1875`
## • `` -> `...1876`
## • `` -> `...1877`
## • `` -> `...1878`
## • `` -> `...1879`
## • `` -> `...1880`
## • `` -> `...1881`
## • `` -> `...1882`
## • `` -> `...1883`
## • `` -> `...1884`
## • `` -> `...1885`
## • `` -> `...1886`
## • `` -> `...1887`
## • `` -> `...1888`
## • `` -> `...1889`
## • `` -> `...1890`
## • `` -> `...1891`
## • `` -> `...1892`
## • `` -> `...1893`
## • `` -> `...1894`
## • `` -> `...1895`
## • `` -> `...1896`
## • `` -> `...1897`
## • `` -> `...1898`
## • `` -> `...1899`
## • `` -> `...1900`
## • `` -> `...1901`
## • `` -> `...1902`
## • `` -> `...1903`
## • `` -> `...1904`
## • `` -> `...1905`
## • `` -> `...1906`
## • `` -> `...1907`
## • `` -> `...1908`
## • `` -> `...1909`
## • `` -> `...1910`
## • `` -> `...1911`
## • `` -> `...1912`
## • `` -> `...1913`
## • `` -> `...1914`
## • `` -> `...1915`
## • `` -> `...1916`
## • `` -> `...1917`
## • `` -> `...1918`
## • `` -> `...1919`
## • `` -> `...1920`
## • `` -> `...1921`
## • `` -> `...1922`
## • `` -> `...1923`
## • `` -> `...1924`
## • `` -> `...1925`
## • `` -> `...1926`
## • `` -> `...1927`
## • `` -> `...1928`
## • `` -> `...1929`
## • `` -> `...1930`
## • `` -> `...1931`
## • `` -> `...1932`
## • `` -> `...1933`
## • `` -> `...1934`
## • `` -> `...1935`
## • `` -> `...1936`
## • `` -> `...1937`
## • `` -> `...1938`
## • `` -> `...1939`
## • `` -> `...1940`
## • `` -> `...1941`
## • `` -> `...1942`
## • `` -> `...1943`
## • `` -> `...1944`
## • `` -> `...1945`
## • `` -> `...1946`
## • `` -> `...1947`
## • `` -> `...1948`
## • `` -> `...1949`
## • `` -> `...1950`
## • `` -> `...1951`
## • `` -> `...1952`
## • `` -> `...1953`
## • `` -> `...1954`
## • `` -> `...1955`
## • `` -> `...1956`
## • `` -> `...1957`
## • `` -> `...1958`
## • `` -> `...1959`
## • `` -> `...1960`
## • `` -> `...1961`
## • `` -> `...1962`
## • `` -> `...1963`
## • `` -> `...1964`
## • `` -> `...1965`
## • `` -> `...1966`
## • `` -> `...1967`
## • `` -> `...1968`
## • `` -> `...1969`
## • `` -> `...1970`
## • `` -> `...1971`
## • `` -> `...1972`
## • `` -> `...1973`
## • `` -> `...1974`
## • `` -> `...1975`
## • `` -> `...1976`
## • `` -> `...1977`
## • `` -> `...1978`
## • `` -> `...1979`
## • `` -> `...1980`
## • `` -> `...1981`
## • `` -> `...1982`
## • `` -> `...1983`
## • `` -> `...1984`
## • `` -> `...1985`
## • `` -> `...1986`
## • `` -> `...1987`
## • `` -> `...1988`
## • `` -> `...1989`
## • `` -> `...1990`
## • `` -> `...1991`
## • `` -> `...1992`
## • `` -> `...1993`
## • `` -> `...1994`
## • `` -> `...1995`
## • `` -> `...1996`
## • `` -> `...1997`
## • `` -> `...1998`
## • `` -> `...1999`
## • `` -> `...2000`
## • `` -> `...2001`
## • `` -> `...2002`
## • `` -> `...2003`
## • `` -> `...2004`
## • `` -> `...2005`
## • `` -> `...2006`
## • `` -> `...2007`
## • `` -> `...2008`
## • `` -> `...2009`
## • `` -> `...2010`
## • `` -> `...2011`
## • `` -> `...2012`
## • `` -> `...2013`
## • `` -> `...2014`
## • `` -> `...2015`
## • `` -> `...2016`
## • `` -> `...2017`
## • `` -> `...2018`
## • `` -> `...2019`
## • `` -> `...2020`
## • `` -> `...2021`
## • `` -> `...2022`
## • `` -> `...2023`
## • `` -> `...2024`
## • `` -> `...2025`
## • `` -> `...2026`
## • `` -> `...2027`
## • `` -> `...2028`
## • `` -> `...2029`
## • `` -> `...2030`
## • `` -> `...2031`
## • `` -> `...2032`
## • `` -> `...2033`
## • `` -> `...2034`
## • `` -> `...2035`
## • `` -> `...2036`
## • `` -> `...2037`
## • `` -> `...2038`
## • `` -> `...2039`
## • `` -> `...2040`
## • `` -> `...2041`
## • `` -> `...2042`
## • `` -> `...2043`
## • `` -> `...2044`
## • `` -> `...2045`
## • `` -> `...2046`
## • `` -> `...2047`
## • `` -> `...2048`
## • `` -> `...2049`
## • `` -> `...2050`
## • `` -> `...2051`
## • `` -> `...2052`
## • `` -> `...2053`
## • `` -> `...2054`
## • `` -> `...2055`
## • `` -> `...2056`
## • `` -> `...2057`
## • `` -> `...2058`
## • `` -> `...2059`
## • `` -> `...2060`
## • `` -> `...2061`
## • `` -> `...2062`
## • `` -> `...2063`
## • `` -> `...2064`
## • `` -> `...2065`
## • `` -> `...2066`
## • `` -> `...2067`
## • `` -> `...2068`
## • `` -> `...2069`
## • `` -> `...2070`
## • `` -> `...2071`
## • `` -> `...2072`
## • `` -> `...2073`
## • `` -> `...2074`
## • `` -> `...2075`
## • `` -> `...2076`
## • `` -> `...2077`
## • `` -> `...2078`
## • `` -> `...2079`
## • `` -> `...2080`
## • `` -> `...2081`
## • `` -> `...2082`
## • `` -> `...2083`
## • `` -> `...2084`
## • `` -> `...2085`
## • `` -> `...2086`
## • `` -> `...2087`
## • `` -> `...2088`
## • `` -> `...2089`
## • `` -> `...2090`
## • `` -> `...2091`
## • `` -> `...2092`
## • `` -> `...2093`
## • `` -> `...2094`
## • `` -> `...2095`
## • `` -> `...2096`
## • `` -> `...2097`
## • `` -> `...2098`
## • `` -> `...2099`
## • `` -> `...2100`
## • `` -> `...2101`
## • `` -> `...2102`
## • `` -> `...2103`
## • `` -> `...2104`
## • `` -> `...2105`
## • `` -> `...2106`
## • `` -> `...2107`
## • `` -> `...2108`
## • `` -> `...2109`
## • `` -> `...2110`
## • `` -> `...2111`
## • `` -> `...2112`
## • `` -> `...2113`
## • `` -> `...2114`
## • `` -> `...2115`
## • `` -> `...2116`
## • `` -> `...2117`
## • `` -> `...2118`
## • `` -> `...2119`
## • `` -> `...2120`
## • `` -> `...2121`
## • `` -> `...2122`
## • `` -> `...2123`
## • `` -> `...2124`
## • `` -> `...2125`
## • `` -> `...2126`
## • `` -> `...2127`
## • `` -> `...2128`
## • `` -> `...2129`
## • `` -> `...2130`
## • `` -> `...2131`
## • `` -> `...2132`
## • `` -> `...2133`
## • `` -> `...2134`
## • `` -> `...2135`
## • `` -> `...2136`
## • `` -> `...2137`
## • `` -> `...2138`
## • `` -> `...2139`
## • `` -> `...2140`
## • `` -> `...2141`
## • `` -> `...2142`
## • `` -> `...2143`
## • `` -> `...2144`
## • `` -> `...2145`
## • `` -> `...2146`
## • `` -> `...2147`
## • `` -> `...2148`
## • `` -> `...2149`
## • `` -> `...2150`
## • `` -> `...2151`
## • `` -> `...2152`
## • `` -> `...2153`
## • `` -> `...2154`
## • `` -> `...2155`
## • `` -> `...2156`
## • `` -> `...2157`
## • `` -> `...2158`
## • `` -> `...2159`
## • `` -> `...2160`
## • `` -> `...2161`
## • `` -> `...2162`
## • `` -> `...2163`
## • `` -> `...2164`
## • `` -> `...2165`
## • `` -> `...2166`
# define and solve prioritization problem
prob <- problem(units, pfeatures) %>%
      add_max_features_objective(budget = prot_area + 10) %>%
      add_feature_weights(weights = ds$tree$edge.length) %>%
      add_relative_targets(0.1) %>%
      add_binary_decisions() %>%
      add_locked_in_constraints(locked_in = protected) %>%
      add_rsymphony_solver(gap = 0)
soln <- solve(prob)
plot(soln + protected, 
     col = c("gray", "orange", "darkgreen"))

Exercises:

  • These exercises were based on the chronogram, which we packaged up with the California flora dataset. Replicate them with the phylogram and cladogram. Do the conservation priorities change?
  • prioritizr has functionality to consider many additional factors, such as land cost, reserve shape, and reserve connectivity. A simple example is add_neighbor_constraints(k = 2), which would force proposed reserves to have at least two neighbors. Like other constraints this will likely add a lot of computational time, so try it when you have a few minutes to wait.
prob <- problem(units, features) %>%
      add_relative_targets(0.1) %>%
      add_binary_decisions() %>%
      add_rsymphony_solver(gap = 0)

soln1 <- prob %>%
      add_max_phylo_div_objective(budget = prot_area,
                                  tree = chronogram) %>%
      solve()
soln2 <- prob %>%
      add_max_phylo_div_objective(budget = prot_area,
                                  tree = cladogram) %>%
      solve()
soln3 <- prob %>%
      add_max_phylo_div_objective(budget = prot_area,
                                  tree = phylogram) %>%
      solve(run_checks = F)

solns <- stack(soln1, soln2, soln3)

# plot priorities for each facet
solns %>%
      setNames(c("chronogram", "cladogram", "phylogram")) %>%
      plot(main = "cladogram",
           col = c("gray", "red"))

# plot number of facets under which each cell is a priority
solns %>%
      sum() %>%
      plot(col = c("gray", "yellow", "orange", "red"))

And that is all!